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We present empirical potentials for dilute transition metal solutes in a-iron. It is in the Finnis- 
Sinclair form and is therefore suitable for billion atom molecular dynamics simulations. First 

S principles calculation shows that there are clear trends across the transition metal series which 
enable us to relate the rescaling parameters to principal quantum number and number of d 
electrons. 

. 

The potential has been developed using a rescaling technique to provide solute-iron and 
O , solute-solute interactions from an existing iron potential. 
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C> 1. Introduction 



The ability to model many component alloys of iron on the atomic level would 
provide an extremely powerful tool for research into the behaviour of these materi- 
als. In particular it would allow the effects of varying proportions of solutes on the 
properties of these materials to be studied in detail. The theoretical insights gained 
would complement the already extensive understanding of such systems found from 
experiment and ultimately impact on the design of new materials to suit particular 
applications such as those for the nuclear industry. 

In principle one would wish to model these materials using ab-initio electronic 
structure calculations but such techniques are prohibitive, being restricted to less 
that 1000 atoms and femto second timescales. Alternative higher level modelling 
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techniques exist, such as kinetic monte-carlo (kMC) and molecular dynamics sim- 
ulations (MD) using empirical potentials, that remove these restrictions but at the 
expense of requiring input from experiment or ab-initio calculations to fix their 
free parameters. In particular, empirical potentials allow billion atom simulations 
to be performed over nano second timescales. The results from such simulations 
are readily used as input to continuum engineering models and ultimately in the 
design of new materials. 

The Finnis-Sinclair scheme [l] was based around the idea of a second moment 
model to the local density of states. In this model the band energy depends on the 
width of the band, the shape of the band, and the occupation of the band. The 
moments theorem 0] shows how band width can be determined from the sum of 
squares of hopping integrals, which forms the physical basis for the cohesive term in 
Finnis-Sinclair potentials. The band shape and occupation are implicitly assumed 
to be constant. 

For elements in a single phase with charge neutrality 0], the d-band occupation is 
essentially constant, and the band shape does not change massively. This under- 
lies the success of single-element potentials. Fitting to alloys has a more troubled 
history. Whereas isoelectronic alloys for isostructural elements work well [J], poten- 
tials for systems involving structural phase transitions or elements from different 
series tend to have poor transferrability from the composition at which they are 
fitted. Thus the generalization to multicomponent alloys seems to require physics 
beyond the second-moment model. Difficulties also arise in the Finnis-Sinclair and 
related schemes such as the embedded atom method (EAM) 0-0] due to the in- 
creasing complexity of the model. For an ^-component system the number of 
parameterised functions grows as N 2 , as does the number of data points required 
to fit the parameters of these functions. 

There is one important exception where the bandshape and electron density are 
reasonably constant, and we might hope that the second-moment approach will 
work. That is multicomponent alloys with one dominant element and multiple 
minority elements, a particularly relevant case here being steel. The dominant 
element fixes both the crystal structure and the electron density and should, in 
principle, connect the behaviour of single solute atoms to that of the dominant 
element and introduce stronger connections between solute-solute interactions and 
the properties of single solutes. 

In this paper we present empirical potentials for single transition metal solutes 
in a-iron by rescaling the functions of a pure iron potential. We also investigate 
ways to connect the interactions between solute particles to those of single solute 
atoms in iron which would allow multi-component potentials to be built once the 
interactions of single solutes in iron are known. We start by providing motivation 
for this procedure from the results of a recent ab-initio study [8[ and then give a 
detailed description of the rescaling strategy. We then discuss the results for single 
solute atoms in iron and present the findings of our investigation into solute-solute 
interactions. Finally we present our conclusions. 
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2. Rescaling 



2.1. Ab-initio calculations: Is rescaling credible? 



If the key physics of substitutional atoms in steel is such that a rescaling approach 
will work, we should expect that the rescaling will involve the d-electron density and 
the principle quantum number. Such an approach should work both for the perfect 
lattice and for defects. The properties of substitutional transition elements in Fe, 
in particular their magnetic character, have long been known to have systematic 
trends 0, [l(J • Here we supplement this work with an emphasis on total energy 
calculations for substitutional atoms and their interactions with point defects in 
bcc Fe. 



We use the VASP code [ll|] with projector augmented wave (PAW) pseudopoten- 
tials [12] and the generalized gradient approximation with the Vosko-Wilk- 
Nusair interpolation lj], which we find to give the best compromise between com- 
putation speed and accuracy. This gives a lattice parameter for pure iron of 2.83A, 
which was used in the impurity and defect calculations to define a fixed- volume su- 
percell. Supercells of 1281 atoms were used with a Monkhorst-Pack 3x3x3 k-point 
grid sampling the Brillouin zone. The energy cutoff was set to 300 eV. Full details 
of the calculations will be published elsewhere [§]. 



The following definition has been used for the binding energy of n defects and 
impurities, {Ai}: 



i=l 



[E({Ai}) + (n - 1)E ] , (1) 



where E(A{) is the energy for a configuration containing Ai only, E{{A{\) refers to a 
configuration containing all the interacting entities and Eq refers to a configuration 
containing no defects or impurities i.e. bulk a-iron. 

These total energy calculations show that there are systematic trends across the 
transition metal series for the free atom substitutional energy, E s , excess pres- 
sure from a single solute, P, first nearest neighbour solute- iron separation, ri nn , 
solute-solute interactions, binding energies of a single solute to a vacancy defect 
at Inn, E^' lnn , and 2nn, E^' 2nn , separations and the binding energies to a (110)- 
self-interstitial defect in the mixed, £? b I,M , compressive, E,', and tensile sites 
^ I,T , (see FIG. [T]for configurations) as shown in FIG. [2j The free atom substitu- 
tional energy was calculated from our ab-initio results for the substitution energy 
from the pure equilibrium phase and the experimental cohesive energies of the pure 
phases |15|] . We take these values as fit targets in order to determine the parameters 
of our potentials, as discussed in the following section. 

We plot all energies against the number of d-electrons in the free atom. In the solid 
this number will be affected by s — d transfer of approximately 0.5 electrons per 
atom. Thus although there are clearly different trends for more-than or less-than 
half filled bands, rigorously defining which material corresponds to a half-filled 
(i-band is not straightforward. 
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(a) (b) 



Figure 1. (a) First and second nearest neighbour solute sites (black) relative to a vacancy defect (square) 
(b) Mixed (M), compressive (C) and tensile (T) solute sites (black) relative to a (HO)-self-interstitial 
defect. 

Two elements produce outlier behaviour: chromium and manganese. Curiously, iron 
chromium and manganese were elements which Finnis Sinclair were unable to fit 
with their original scheme [lj . These elements exhibit unusual magnetic behaviour, 
which presumably accounts for this. 



2.2. Rescaling Strategy 



The starting point for our fitting strategy is the pure iron EAM potential of Ackland 
et al. [la]. We have chosen this iron potential over those from more recent works 0- 
[lj| because it reproduces many of the properties of iron despite its relatively simple 
form. 

The most general form for the energy, U, of an EAM potential is given by 

U({r ab }) = V (Xa ' Xb) (r ab ) + ^F^(p a ), (2) 

a 

(3) 



a,b>a 



where V^ Xa,Xh \ (f)( Xa > Xb > and F^ Xa ^ are parameterised functions dependent on the 
atomic species, X a and X b . The cross-species pair functions are taken to be sym- 
metrical here, i.e. V^ X ' Y ^ = V^ Y,X ' when X 7^ Y, as are the functions, ^ X,Y \ 

We use the same forms for the component functions of our potential as used in the 
pure iron potential [lf| • In particular we define the pair functions by 



V (X,Y) {7 



47reo'' 

exp(S + Bir + B 2 r 2 + B 3 r 3 ) 



lc?(W(r) = Ej=i 



(X,Y),{X,Y) 



r) 3 H(r 



(X,Y) 



r <r\ 
r) r > T2 



(4) 

where Zx is the atomic number of species X, r s = 0.88534 a b /JZ 2 J 3 + Zy /3 , db IS 
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Figure 2. Fit targets (squares and lines) and corresponding values from our empirical model (circles) versus 
number of d-electrons for 3d solutes (black), 4d solutes (red) and 5d solutes (green). The first nearest 
neighbour solute-iron separation, f^~ Fe , is plotted relative to the corresponding iron-iron separation, 

Fe-Fe 
' Inn 
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the Bohr radius and 

£(x) = 0.1818e- 3 - 2a; + 0.5099e- a9423a; + 0.2802 e -°- 4029a; + 0.02817e-°- 2016x . (5) 



The functional form used below n = 0.9A is the universal screened potential 
of Biersack and Ziegler [13], above r 2 = 1.9 A is a parameterised cubic spline with 
cutoffs implemented by the use of Heaviside step functions, H, and between these is 
an interpolating function that ensures continuity of the function and its derivative. 

For the embedding functions, F^ x \ we take the standard square root form for all 
atomic species i.e. F^ x \p) = —Jp~. 

The cj> functions take the form of a simple cubic spline 

^\r) = Y.A^ Y \R^ - r?H(R^ - r). (6) 
k=l 



For the pure iron component functions i.e. y( Fe ' Fe ) and 0( Fe > Fe ) we take the param- 
eters directly from Iron-solute interactions are defined by rescaling these two 
functions using rescale parameters, {pi }: 

C^ x \p[ X) r)=p { 2 X) C^ Fc \r) (7) 

(Fe,X) (p (X) r)=p (X) 0(Fe , F e) (r) (g) 

This is equivalent to a direct rescaling of the parameters of the cubic spline func- 
tions given by, for example, 



,(Fe,X) _ (X) (Fe,Fe) 



(9) 



<r x > - ^L"f- re) - (id 

(X) 

We take the rescaling factors, {p\ }, to be the adjustable parameters for the 
purposes of fitting. The trends in the fit target data should therefore translate 
to trends in these rescale parameters across the transition metal series. In fact it 
should be possible to quantify these trends by finding functional forms that relate 
the rescale parameters to the elementary electronic properties of the solutes. We 

present our results for the rescale parameters and a functional form for them in 

(x) 

terms of the number of d-electrons per atom, n\ , in the following section. 
Solute-solute interactions are defined by a similar rescaling procedure: 

C^ ¥ \p[ X ' Y) r) = p { 2 X ' Y) C^ Fc \r) (11) 
^ x ' Y \p (X ' Y) r) = p[ X ' r V (Fe ' Fe) (r)- (12) 



However, we do not determine these rescale parameters from fitting. Instead we 
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relate them to the rescale parameters for the iron-solute interactions, i.e. 

(x,Y) (x,Y)r PO 0O\ f-io\ 
Pi =Pl (Pi ,Pi )• (13) 

This is the key step that ensures we can construct multi-component alloys once 
the iron-solute interactions are known. 



3. Single solute interactions in iron 



(X) 

In order to determine the rescale parameters, {p\ }, we fit to the ab-initio data 
shown in FIG. [2] for transition metal solutes in a-iron. We do not, however, fit to 
E^ l ' M as no satisfactory results were found with this quantity included. 

The fitting procedure was accomplished by minimising a standard least squares 
response function, % 2 , of the fit parameters, {pi}, given in terms of the fit targets, 
{t r }, model values, {m r ({pi})}, and weight factors, {cr r }, by 

Our potential model values were all calculated in atomically relaxed 4x4x4 bcc unit 
cell configurations, i.e. 128 atoms before the introduction of defects and solutes, 



at the equilibrium volume for the pure iron potential, i.e. ao = 2.8665A 16]. This 
was done in order to appropriately match the ab-initio fit target data. We chose 
weight factors of O.OleV (or 1% of the fit target value if that is larger) for energies, 

0.005 A for lengths and 5 x 10~ 4 eV/A for pressures in our fits. 

The fitted rescale parameters are given in FIG.[3]and TABLE. IA1I It it immediately 
clear that there are trends across the series, especially for the 4d and 5d transition 
metal solutes. The one notable exception is manganese whose anomalous rescale 
parameters match the equally anomalous properties of the solute itself. 

Looking at FIG. [2] we can see that the trends in the rescale parameters translate 
to the potential model values themselves. There is especially good reproduction 
of the substitution energy across all three series. The binding energy, E^' , is 
reproduced almost as well but shows slight deviation from the ab-initio target data 
at the ends of the series and especially for low rid- Such deviations from the usually 
parabolic trends in the ab-initio data are seen ^generally for the other fit targets. 
The most notable is for the binding energy, E, 1 ' , where the potential models show 
approximately linear behaviour in n<j. The overall deviations and the quality of the 
fits is best quantified via the response function values, as shown in TABLE. IA2I It 
is clear from this data that the potentials for the low rid solutes perform especially 
poorly. This is an interesting results because the influence of s-electrons become 
increasingly important for these elements and their effects are not included in the 
pure iron potential we have rescaled here. 



It is also clear from the response function data that the 3d solute potentials perform 
better overall than those of the other two series despite the presence of complex 
magnetic interactions. Even the anomalous properties of manganese are repro- 
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Figure 3. Rescale parameters as a function of the number of d-electrons for the 3d solutes (black circles), 

(X) (X) (X) 

4d solutes (red squares) and 5d solutes (green diamonds). The graphs are for (a) , (b) p 2 , (c) p^ 
and (d) p. . 

duced well. The excess solute pressures are, however, significantly underestimated 
although this is true of the 4d and 5d elements also. 



Finally it is worth returning to the difficulties experienced in fitting the binding 



energy of the mixed interstitial, E b ' . As can be seen from FIG. 2(f) our potentials 
significantly underestimate the magnitude of this value for the 4d and 5d elements. 
Including this value in the fits did result in a more accurate reproduction but at 
too much cost to the reproduction of the other fit targets. Despite this the mixed 
site is still preserved as the least favoured for a solute to occupy around a (110)- 
self-interstitial defect. This failure is unlikely to have serious consequences for any 
molecular dynamics simulations using our potentials, since the mixed dumbbell site 
is not a migration barrier. The interstitial is repelled by the solute by well above 
any realistic thermal energy and so the mixed interstitial site will occur with very 
low probability. 



4. Solute-solute interactions in iron 



In order to gain some insight into the possible relationship between the solute-solute 
rescale parameters and the solute-iron rescale parameters we have fitted the rescale 
parameters, {p\ X ' X ^}, to reproduce ab-initio values [H] for the solute-solute binding 
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(c) (d) 

Figure 4. Rescale parameters as a function of the number of d-electrons for the 3d solutes (black circles), 

( X X ) ( X X ) 

4d solutes (red squares) and 5d solutes (green diamonds). The graphs are for (a) p\ ' , (b) pi, ' > ( c ) 
Pg ' and (d) p\ . 

energies from Inn to 5nn separation, EL - ,inn , and to the separations between 
solutes at Inn separation, and 2nn separation, The resulting values 

for the rescale parameters are shown in FIG. [Hand our model values are compared 
with the ab-initio fit targets in FIG. [SJ 

The fit targets in FIG. [5] once again show a two-part trend across the group, which 
gives us hope that not only can rescaling be used to fit the potential functions di- 
rectly, but also that the rescale parameters can be deduced directly from the atomic 
number. The picture emerging from the fit parameters themselves FIG. [4] is less 
clear. For elements above half-filling there are clear trends with principle quantum 
number and rid- However, for 5d elements with rid < 5 there is considerable scatter 
for p% X ^ and p^ ■ It appears that for a longer-ranged <j) can compensate for a 
stronger repulsion, since this anomaly is not present in the fit targets in FIG. [5J it 
must be an artifact of the fitting process itself. 



5. Conclusions 



In conclusion, we have advanced the hypothesis that the interactions between tran- 
sition metal atoms in iron can be described by simply scaling the parameters of a 
Finnis-Sinclair model, and further that these scaling parameters are simply func- 
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(c) (d) 

Figure 5. Solute-solute fit targets (squares and lines) and corresponding values from our empirical model 
(circles) versus number of d-electrons for 3d solutes (black), 4d solutes (red) and 5d solutes (green). 

tions of the number of (i-electrons and the principal quantum number. Quantum 
mechanical calculations of interactions between solutes and point defect show clear 
trends across the group. We have presented best-fit rescaled potentials for all tran- 
sition metal elements in an iron matrix. 

Our hypothesis is based on the notion that the environment around the impurity 
is close to that in magnetic bcc iron, hence we do not expect that the rescaled 
potentials will be transferrable to very different electronic environments such as 
pure elements. 
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Solute, X 



(X) 

Pi 



(X) 

P2 



(X) 
P3 



(X) 

Pi 
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Tabic A2. Response function values, x 2 , for the best rescale parameters given in TABLE. IaTI 



Solute,X 


X 2 (best) 


Ti 


568.9656946430212 


V 


221.2048362240138 


Cr 


99.63599215304797 


Mn 


178.81158941110317 


Co 


43.901666493113794 


Ni 


65.72761363633806 


Cu 


204.8091878556158 


Zr 


4796.903769001561 


Nb 


1689.7527370510506 


Mo 


216.13918658727746 


Tc 


137.2039869893378 


Ru 


69.61819599245268 


Rh 


26.410945708888832 


Pd 


103.65177965112031 


Ag 


477.68739615717334 


Hf 


3491.8432724139166 


Ta 


1502.974429434891 


W 


408.3356703154182 




1 98 776969870386 


Os 


85.53222566998141 


Ir 


35.81070622123357 


Pt 


69.88189007267775 


Au 


386.0176534820233 



